home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / U-Z / VideoToolBox Folder / VideoToolboxSources / Luminance.c < prev    next >
Encoding:
Text File  |  1993-03-04  |  36.9 KB  |  991 lines  |  [TEXT/KAHL]
  1. /*
  2. Luminance.c
  3. Version 1.95
  4. See Luminance.h for prototypes.
  5. See Luminance.doc for explanation.
  6. Also see:
  7. D.G. Pelli and L. Zhang (1991) Accurate control of contrast on microcomputer displays. 
  8. Vision Research, 31:1337-1360.
  9.  
  10. Copyright (c) 1989-1993 Denis G. Pelli. This software is free; you may use it in
  11. your research and give it away to others, with the following restrictions. Any
  12. copy you give away must include this paragraph, unmodified, and, if you have
  13. made any changes to the rest of the file, must include a note, added to HISTORY
  14. below, giving your name, the date, and a description of the changes. This
  15. software may not be sold, whether in source or compiled form, without my
  16. permission. I hope you will find this software useful, but I can't promise that
  17. it will work for you, and am not offering any support. That's why it's free. I
  18. would appreciate reports of bugs and improvements.
  19.  
  20. Denis G. Pelli, Ph.D.
  21. Institute for Sensory Research, Syracuse University, Syracuse, NY 13244-5290, U.S.A.
  22. denis_pelli@isr.syr.edu
  23.  
  24. HISTORY:
  25.  
  26. 4/24/89 dgp first working version.
  27. 4/28/89 dgp added device argument.
  28. 6/23/89 dgp added support for a video attenuator that has unequal gain in the three
  29.             channels. Added LoadClut() and LuminanceClut() routines.
  30. 7/30/89 dgp replaced call to SetEntries (which calls the Color Manager) by a call
  31.             to GDSetEntries (which calls the device driver).
  32. 8/10/89 dgp fixed bugs.
  33. 8/13/89 dgp rewrote LToV() to use nth order polynomial. Broke up LinearizeClut into
  34.             three routines--SetLuminanceRange, SetLuminance, and SetLuminances--the
  35.             last of which is equivalent to the old LinearizeClut.
  36. 8/16/89    dgp    After reading the Brooktree DAC manual I changed the documentation and
  37.             increased tolerance from 0.5 to 1 LSB, and now use the highest
  38.             luminance gain (worst case) over the requested range to transform it
  39.             into a luminance difference tolerance.
  40. 8/21/89    dgp    Made minor improvements in LToV() to deal with failure to converge due to 
  41.             getting stuck in a local minimum of the quartic. My solution is to now
  42.             REQUIRE that Lmin and Lmax be filled in by the user in the calibration
  43.             structure. In practice this will guarantee that the minimum luminance
  44.             will not be below the local minimum found at the base of the rising curve.
  45.             Last week I also changed the LToV algorithm slightly, to do a bisection if
  46.             Newton's method would take us out of range.
  47. 9/7/89 dgp    Fixed bug in sorting routine, Sort().
  48. 9/10/89 dgp    Commented out the polynomial versions of VToL and LToV and replace them by
  49.             new versions that use a power law. The power law is a marginally better fit
  50.             than a 6th order polynomial (using only 4 instead of 7 parameters) and its
  51.             inverse can be computed much more quickly, about 0.3 ms instead of 2 ms.
  52. 9/11/89 dgp    Deciding to be indecisive, I introduced a conditional POWER_LAW_FIT
  53.             to control the choice of power law or polynomial fit for LToV and VToL.
  54. 9/16/89    dgp    Having finished writing the paper (Pelli & Zhang, 1991) I came back to
  55.             change the program to agree with what I wrote. I changed the equivalent
  56.             number tolerance to be the sum instead of the maximum of the 
  57.             variable-dac gains.
  58. 9/25/89 dgp Fixed minor bug in computation of tolerance. Now gives correct answer even
  59.             when the lowLuminance or highLuminance is out of range.
  60. 10/9/89 dgp    Made minor change so that when some of the gain[] are zero SetLuminance will
  61.             load zero into the unused colorSpec array elements. This is only a cosmetic
  62.             change.
  63. 10/28/89 dgp Made minor changes. I declared LToVPolynomial and LToVPower and
  64.             VToLPolynomial and VToLPower instead of having a conditional compilation. 
  65.             This makes both flavors of routine always available. I also changed 
  66.             LToVPolynomial to use LToVPower instead of LToVQuadratic for its initial 
  67.             guess. The quadratic guess was often poor and occasionally led to 
  68.             convergence problems. If the power law fit is not available then it reverts
  69.             to using the quadratic fit. If that's not available then it reverts to using
  70.             a middle-of-range guess.
  71. 10/28/89 dgp I eliminated the compile-time flag POWER_LAW_FIT and instead use whichever, 
  72.             polynomial or power, provides a better fit, as determined at run time. I 
  73.             biased the comparison of rms error to favor the power law fit since it has 
  74.             few parameters and can be inverted much more quickly. Note: if the 
  75.             power, polynomial, or quadratic fit parameters are not supplied then the 
  76.             appropriate LR.powerError, LR.polynomialError, or LR.quadraticError field 
  77.             should be set to infinity: INF or 1.0/0.0.
  78. 12/5/89    dgp    Added the trivial routine LToL() which enforces the bounds LP->LMin and 
  79.             LP->LMax.
  80. 4/2/90    dgp    Capitalized the word "To" in all function names, e.g. LToL().
  81. 4/22/90    dgp    Version 1.4. Made minor changes to the documentation. Renamed LToV() to
  82.             LToVFormulaic(), and introduced a new LToV() that uses a table to run faster.
  83.             LToV() takes an average of 170 µs, whereas LToVFormulaic() takes
  84.             1700 µs. I also now use the luminanceTable, if available, to speed up
  85.             VToL() slightly, from 310 µs to 180 µs. Each call to SetLuminance() now
  86.             takes 0.9 ms whereas it used to take 2.5 ms. The loss in accuracy is
  87.             negligible. The interpolation error can be reduced still further by
  88.             increasing LUMINANCES_IN_TABLE. Each doubling of the table size quarters the
  89.             maximum possible error of the interpolation.
  90. 4/23/90    dgp    Added a few lines of code to LToV() to check near the last index, in the
  91.             spirit of Numerical Recipes in C, hunt.c, before starting the bisection
  92.             search. I updated the times in the 4/22/90 note to reflect the latest timing
  93.             by TestLuminance.
  94. 7/27/90    dgp    Tightened up the error checking for illegal entries into the clut. I have
  95.             the impression that, contrary to specification, the Apple video driver
  96.             crashes and corrupts itself if given an out-of-range entry value in a
  97.             setEntry Control call.
  98. 9/18/90    dgp    Changed all instances of "v" to "V". The final version of the Pelli
  99.             & Zhang (1991) paper refers to the nominal voltage v; this file
  100.             now refers to the "equivalent number" V; they are related by V=255*v.
  101. 9/24/90    dgp    Updated the documentation to correspond more closely to the (hopefully)
  102.             final version of the Pelli & Zhang (1991) paper.
  103. 10/29/90 dgp Doubled speed of SetLuminance(). Now SetLuminance() takes 133 ms to
  104.             do a complete clut for a small new luminance range, and takes 188 ms for a
  105.             large new range. This required minor changes to SetLuminanceRange() & LToV().
  106.             Changed all function headers to ANSI prototype style.
  107. 10/30/90 dgp Replaced phrase "clut value" by "nominal voltage", for consistency with
  108.             Pelli & Zhang (1991).
  109.             Introduced new fixed fraction data type called "Milli", defined in Milli.h.
  110.             By judicious replacement of double by Milli, particularly in the luminance
  111.             table, I have increased the speed of SetLuminance by a factor of three.
  112.             A Mac II now takes only 42 to 50 ms to build a whole clut. This new feature
  113.             is enabled by setting FAST_LUMINANCE to 1 in Luminance.h. There is a slight
  114.             increase in the luminance tolerance.
  115. 10/31/90 dgp More fine tuning. Now takes 31 to 44 ms to build a whole clut. Introduced
  116.             conditional FAST_LUMINANCE in Luminance.h that causes the same code to
  117.             be compiled with either Milli or double variables. The _SetLuminance()
  118.             routine is now quite polished, and runs fast. All variables and routines
  119.             whose names begin with underscore _ are written here with type Milli, but
  120.             if FAST_LUMINANCE is false then this type is redefined (in this file only)
  121.             as double, and all the Milli macros are redefined to be appropriate for
  122.             a double. So bear in mind that things beginning with underscore are of
  123.             unknown type. The virtue of being able to run the same code with Milli or
  124.             double computations is that if you suspect an error in the (homemade)
  125.             Milli arithmetic, you can try out double arithmetic simply by setting
  126.             FAST_LUMINANCE to 0 and recompiling.
  127.             In the interest of speed I have streamlined the tolerance
  128.             computation. The new tolerances are probably as good as the old ones.
  129. 11/2/90 dgp Added _SetLuminances() subroutine that substitutes linear interpolation
  130.             for most of the calls to _LToV(). This has greatly speeded up
  131.             SetLuminancesAndRange(), with minor loss of accuracy. A complete clut now
  132.             takes 9 to 31 ms. Even a Mac II can now compute low-contrast cluts
  133.             on the fly, between frames. (The user may wish to optimize the speed-
  134.             accuracy trade-off controlled by the LINEAR_V_DOMAIN parameter, which
  135.             determines the maximum width of the interpolation interval.)
  136. 11/6/90 dgp Replaced Milli by FIXED, which can be compiled as either double or Fixed.
  137.             The Fixed math is sufficiently precise as to give the same DAC values as
  138.             double math.
  139.             Now figure out optimal bit shift LT->L.LShift to preserve as much resolution
  140.             as possible when interpolating in _LToV(). 
  141.             A complete clut now takes 8 to 30 ms. For threshold stimuli the average
  142.             time will usually be less than 15 ms.
  143. 11/8/90 dgp    Eliminated gamma slope table since the speed-up it offered was too small
  144.             to measure. Tidied up the computation of LShift.
  145. 8/24/91    dgp    Made compatible with THINK C 5.0.
  146. 12/17/91 dgp Added new routines: IncrementLuminance() and GetV().
  147.             The former is used to obtain the lowest possible contrast. On a log contrast
  148.             scale, minus infinity (zero contrast) is a poor approximation to even
  149.             a small log contrast.
  150. 3/11/92    dgp    Minor editing of comments.
  151. 12/2/92 dgp Fixed minor bug in SetLuminanceRange() reported by Wei Xie and Ken Alexander
  152.             that could cause the THINK C Debugger to report an error when an 
  153.             uninitialized double V[] was converted to an int. The value was not used 
  154.             for anything, so this was innocuous.
  155. 12/17/92 dgp 1.94 Began enhancement to allow any dacSize, not just 8-bit dacs. I have
  156.             confirmed that everything still works fine with 8-bit dacs, but I don't
  157.             have a 9-bit dac video card to test it any further. •CAUTION: this probably
  158.             does not yet work correctly for 9-bit dacs. E.g. it seems likely that the 
  159.             fixed point math will break if VMax is increased to 511; a lot of careful 
  160.             error analysis went into the original coding, at which time I thought that 
  161.             it was safe to assume an 8-bit dac (i.e. VMax==255). •Removed obsolete 
  162.             support for THINK C 4. 
  163. 12/21/92 dgp No longer load unused dac bits.
  164. 1/4/93    dgp    1.95 LoadLuminances now checks the gdType.
  165. */
  166. #include "VideoToolbox.h"        /* prototype for GDSetEntries() */
  167. #include "Luminance.h"
  168. #include <assert.h>
  169. #include <math.h>
  170.  
  171. static void Sort(int n,double arr[],int krr[]);    /* for internal use only */
  172.  
  173. #undef LongToFix
  174. #undef FixToLong
  175. #undef DoubleToFix
  176. #undef FixToDouble
  177. #if FAST_LUMINANCE
  178.     #define FIXED Fixed
  179.     #if !MACINTOSH
  180.         typedef long Fixed;
  181.         #define FixMul(x,y) DoubleToFix(FixToDouble(x)*FixToDouble(y))
  182.         #define FixDiv(x,y) DoubleToFix(FixToDouble(x)/FixToDouble(y))
  183.     #endif
  184.     #define LongToFix(x) ((long)(x)<<16)
  185.     #define FixToLong(x) ((x)>>16)
  186.     #define DoubleToFix(x) ((Fixed)((x)*65536.+0.5))
  187.     #define FixToDouble(x) ((double)(x)*(1./65536.))
  188.     #define D(x) FixToDouble(x)                            /* handy for debugging */
  189. #else
  190.     #define FIXED double
  191.     #define FixMul(x,y) ((double)(x)*(y))
  192.     #define FixDiv(x,y) ((double)(x)/(y))
  193.     #define LongToFix(x) ((double)(x))
  194.     #define FixToLong(x) ((long)(x))
  195.     #define DoubleToFix(x) ((double)(x))
  196.     #define FixToDouble(x) ((double)(x))
  197. #endif
  198.  
  199. #undef MAX
  200. #undef MIN
  201. #define MAX(a,b) ((a) > (b) ? (a) : (b))
  202. #define MIN(a,b) ((a) < (b) ? (a) : (b))
  203.  
  204. double SetLuminance(GDHandle device,luminanceRecord *LP
  205.     ,int theEntry,double luminance
  206.     ,double lowLuminance,double highLuminance)
  207. /*
  208. Set one entry in the ColorSpec table (and the clut if device is not NULL) to
  209. a specified luminance. It's ok for lowLuminance to be greater than highLuminance; 
  210. they still designate a range.
  211. */
  212. {
  213.     FIXED _luminance;
  214.     
  215.     if(theEntry<0 || theEntry>=COLORS 
  216.         || device!=NULL && theEntry>=GDCOLORS(device))
  217.             PrintfExit("\007SetLuminance: illegal entry %d\n",theEntry);
  218.     
  219.     SetLuminanceRange(LP,lowLuminance,highLuminance);
  220.     _luminance=DoubleToFix(luminance);
  221.     _SetLuminance(LP,theEntry,_luminance);
  222.     if(device != NULL)LoadLuminances(device,LP,theEntry,theEntry);
  223.     return FixToDouble(_Tolerance(LP,_luminance));
  224. }
  225.  
  226. double SetLuminances(GDHandle device,luminanceRecord *LP
  227.     ,int firstEntry,int lastEntry
  228.     ,double firstLuminance,double lastLuminance)
  229. /*
  230. Set a series of entries in the ColorSpec table (and the clut if device is not NULL)
  231. to a linear sequence of luminances. Assume this is the entire luminance range of 
  232. interest.
  233. */
  234. {
  235.     return SetLuminancesAndRange(device,LP,firstEntry,lastEntry
  236.         ,firstLuminance,lastLuminance,firstLuminance,lastLuminance);
  237. }
  238.  
  239. double SetLuminancesAndRange(GDHandle device,luminanceRecord *LP
  240.     ,int firstEntry,int lastEntry
  241.     ,double firstLuminance,double lastLuminance
  242.     ,double lowLuminance,double highLuminance)
  243. /*
  244. Set a series of entries in the ColorSpec table (and the clut if device is not
  245. NULL) to a linear sequence of luminances. Uses last two arguments to set the
  246. luminance range.
  247.  
  248. I introduced a stack-space check before calling _SetLuminances(), since it calls
  249. itself recursively and may use 1000 bytes all together. However, printing an
  250. error message requires at least 4500 bytes. So I quit if the stack gets below
  251. 6000, so that the user will be presented with an understandable error message
  252. rather than a mysterious quit or crash.
  253. */
  254. {
  255.     FIXED _tolerance,_t,_dL64;
  256.     FIXED _firstL,_lastL,_firstV,_lastV;
  257.     
  258.     if(firstEntry<0 || firstEntry>lastEntry || lastEntry>=COLORS 
  259.         || device!=NULL && lastEntry>=GDCOLORS(device) )
  260.         PrintfExit("\007SetLuminancesAndRange: illegal entries %d %d\n",firstEntry,lastEntry);
  261.     if(StackSpace()<6000)PrintfExit("SetLuminancesAndRange: not enough stack space.\007\n");
  262.     SetLuminanceRange(LP,lowLuminance,highLuminance);
  263.     _firstL=DoubleToFix(firstLuminance);
  264.     _lastL=DoubleToFix(lastLuminance);
  265.     _firstV=_LToV(LP,_firstL + LP->_LOffset);
  266.     _lastV=_LToV(LP,_lastL + LP->_LOffset);
  267.     if(lastEntry != firstEntry)
  268.         #if FAST_LUMINANCE
  269.         _dL64=DoubleToFix(64.0*(lastLuminance-firstLuminance)/(lastEntry-firstEntry));
  270.         #else
  271.         _dL64=DoubleToFix((lastLuminance-firstLuminance)/(lastEntry-firstEntry));
  272.         #endif
  273.     else
  274.         _dL64=0;
  275.     _SetLuminances(LP,firstEntry,lastEntry,_firstL,_dL64,_firstV,_lastV);
  276.     if(device != NULL)LoadLuminances(device,LP,firstEntry,lastEntry);    
  277.     
  278.     /* quick and dirty estimate of tolerance */
  279.     _tolerance=LP->L._L[0] - _firstL;
  280.     _t=_lastL - LP->L._L[LUMINANCES_IN_TABLE-1];
  281.     if(_tolerance<_t)_tolerance=_t;
  282.     if(_tolerance<0)_tolerance=0;
  283.     _tolerance+=LP->_tolerance;
  284.     return FixToDouble(_tolerance);    /* estimate of max error in ΔL */
  285. }
  286.  
  287. void _SetLuminances(luminanceRecord *LPtr,int first,int last
  288.     ,FIXED _firstL,FIXED _dL64,FIXED _firstV,FIXED _lastV)
  289. /*
  290. This routine eliminates most of the the slow _LToV() table lookups and
  291. interpolations by assuming that the gamma function is smooth enough that we may
  292. linearly interpolate over any interval of V no wider than LINEAR_V_DOMAIN. Since
  293. the gamma function is roughly parabolic, the consequent error in L will be
  294. proportional to the square of LINEAR_V_DOMAIN, and independent of where this
  295. interval is along V. Since we're inscribing straight line segments in a
  296. positively accelerating gamma function, we will overestimate luminance, and
  297. consequently V will be too low. The error will be greatest at the middle of each
  298. linearly interpolated V interval.
  299.  
  300. If the requested V interval is larger than LINEAR_V_DOMAIN, then it is split
  301. into two intervals by making a pair of recursive calls. LINEAR_V_DOMAIN is
  302. defined in Luminance.h. I suggest a value of 4, but it may be set larger to
  303. attain slightly higher speed at lower accuracy.
  304.  
  305. CAUTION: This routine is a stack hog. Each call uses up about 64 bytes on the
  306. stack, and it calls itself recursively, up to log2(256/LINEAR_V_DOMAIN) times.
  307. Do NOT change any of the declarations of the "new" variables used in the first
  308. if{} to "static", as that would cause the recursion to fail.
  309. */
  310. {
  311.     static int doLast=1;
  312.     int saveDoLast;
  313.     int newOne;
  314.     FIXED _newL,_newV;
  315.     register luminanceRecord *LP=LPtr;
  316.     register int i,Vi;
  317.     register FIXED _VToGo,_g;
  318.     static RGBColor *RGBPtr;    /* static in order to minimize stack usage */
  319.     static int theEntry;        /* static in order to minimize stack usage */
  320.     static FIXED _V,_dV;        /* static in order to minimize stack usage */
  321.     register short leftShift;
  322.     
  323.     if(last-first>1 ){
  324.         _dV=_lastV-_firstV;
  325.         if(_dV>LongToFix(LINEAR_V_DOMAIN) || _dV<LongToFix(-LINEAR_V_DOMAIN)){
  326.             newOne=(first+last)>>1;
  327.             #if FAST_LUMINANCE
  328.                 _newL=_firstL+((newOne-first)*_dL64>>6);
  329.             #else
  330.                 _newL=_firstL+(newOne-first)*_dL64;
  331.             #endif
  332.             _newV=_LToV(LP,_newL + LP->_LOffset);
  333.             saveDoLast=doLast;
  334.             doLast=0;
  335.             _SetLuminances(LP,first,newOne,_firstL,_dL64,_firstV,_newV);
  336.             doLast=saveDoLast;
  337.             _SetLuminances(LP,newOne,last,_newL,_dL64,_newV,_lastV);
  338.             return;
  339.         }
  340.     }
  341.     if(last!=first)_dV=(_lastV-_firstV)/(last-first);
  342.     else _dV=0;
  343.     if(!doLast)last--;
  344.     leftShift=LP->leftShift;
  345.     for(theEntry=first,_V=_firstV;theEntry<=last;theEntry++,_V+=_dV){
  346.         /****** This section of code is copied from _SetLuminance() ****/
  347.         RGBPtr = &LP->table[theEntry].rgb;
  348.         *RGBPtr=LP->rgb;                        /* load for fixed DACs */
  349.         _VToGo=_V - LP->_VFixed + LP->_VHalfStep;
  350.         for(i=LP->fixed;i<LP->dacs;i++) {
  351.             _g=LP->_gain[i];
  352.             Vi=_VToGo/_g;                        /* truncate to integer */
  353.             if(Vi>LP->VMax)Vi=LP->VMax;
  354.             if(Vi<LP->VMin)Vi=LP->VMin;
  355.             _VToGo -= Vi*_g;
  356.             ((short *)RGBPtr)[LP->dac[i]]=Vi<<leftShift;
  357.         }
  358.         /***************************************************************/
  359.     }
  360.     return;
  361. }
  362.  
  363. void _SetLuminance(luminanceRecord *LPtr,int theEntry,FIXED _luminance)
  364. /*
  365. Set one entry in the ColorSpec table to a specified luminance. This is the
  366. private subroutine that actually does all the work for SetLuminance(). This
  367. routine is designed to run as fast as possible.
  368. */
  369. {
  370.     register luminanceRecord *LP=LPtr;
  371.     register int i,Vi;
  372.     register FIXED _VToGo,_g;
  373.     RGBColor *RGBPtr;
  374.     register short leftShift=LP->leftShift;
  375.     
  376.     RGBPtr = &LP->table[theEntry].rgb;
  377.     *RGBPtr=LP->rgb;                        /* load for fixed DACs */
  378.     _VToGo=_LToV(LP,_luminance + LP->_LOffset) - LP->_VFixed + LP->_VHalfStep;
  379.     for(i=LP->fixed;i<LP->dacs;i++) {
  380.         _g=LP->_gain[i];
  381.         Vi=_VToGo/_g;                        /* truncate to integer */
  382.         if(Vi>LP->VMax)Vi=LP->VMax;
  383.         if(Vi<LP->VMin)Vi=LP->VMin;
  384.         _VToGo -= Vi*_g;
  385.         ((short *)RGBPtr)[LP->dac[i]]=Vi<<leftShift;
  386.     }
  387. }
  388.  
  389. void IncrementLuminance(GDHandle device,luminanceRecord *LPtr,int theEntry)
  390. /*
  391. Make smallest possible increase of the luminance of one entry in the ColorSpec
  392. table.
  393. */
  394. {
  395.     double V;
  396.     register luminanceRecord *LP=LPtr;
  397.     
  398.     V=GetV(device,LP,theEntry);
  399.     V+=LP->VHalfStep;
  400.     V+=LP->VHalfStep;
  401.     SetLuminance(device,LP,theEntry,VToL(LP,V),LP->lowLuminance,LP->highLuminance);
  402. }
  403.  
  404. FIXED _Tolerance(luminanceRecord *LPtr,FIXED _luminance)
  405. {
  406.     register luminanceRecord *LP=LPtr;
  407.     register FIXED _tolerance,_t;
  408.     
  409.     /* quick and dirty estimate of tolerance */
  410.     if(LP->L.exists==luminanceSet){
  411.         _tolerance=LP->L._L[0] - _luminance;
  412.         _t=_luminance - LP->L._L[LUMINANCES_IN_TABLE-1];
  413.         if(_tolerance<_t)_tolerance=_t;
  414.         if(_tolerance<0)_tolerance=0;
  415.         _tolerance+=LP->_tolerance;
  416.     }
  417.     else _tolerance=LP->_tolerance;
  418.     return _tolerance;
  419. }
  420.  
  421. double GetLuminance(GDHandle device,luminanceRecord *LP,int theEntry)
  422. /*
  423. If device is not NULL then examines one entry in the actual clut, otherwise
  424. examines the ColorSpec table contained in *LP, and in either case returns the
  425. luminance that will be produced. Supplying an illegal entry value results in a
  426. returned value of -INF.
  427. */
  428. {
  429.     return VToL(LP,GetV(device,LP,theEntry));
  430. }
  431.  
  432. double GetV(GDHandle device,luminanceRecord *LP,int theEntry)
  433. /*
  434. If device is not NULL then examines one entry in the actual clut, otherwise
  435. examines the ColorSpec table contained in *LP, and in either case returns the V
  436. that will be produced. Supplying an illegal entry value results in a returned
  437. value of -INF.
  438. */
  439. {
  440.     RGBColor *myRGBPtr=NULL;
  441.     ColorSpec myCSpec;
  442.     double V;
  443.     int error;
  444.     
  445.     if(device != NULL) {
  446.         if(theEntry<0 || theEntry>=GDCOLORS(device)){
  447.             PrintfExit("GetLuminance: illegal entry %d\n",theEntry);
  448.         }
  449.         error=GDGetEntries(device,theEntry,0,&myCSpec);
  450.         if(error){
  451.             PrintfExit("GetLuminance: GDGetEntries error %d\007",error);
  452.         }
  453.         myRGBPtr=&myCSpec.rgb;
  454.     }
  455.     else {
  456.         if(LP!=NULL && theEntry>=0 && theEntry<COLORS) 
  457.             myRGBPtr=&LP->table[theEntry].rgb;
  458.     }
  459.     if(myRGBPtr == NULL) return -INF;
  460.     V = LP->r*(myRGBPtr->red>>8);
  461.     V += LP->g*(myRGBPtr->green>>8);
  462.     V += LP->b*(myRGBPtr->blue>>8);
  463.     return V;
  464. }
  465.  
  466. void LoadLuminances(GDHandle device, luminanceRecord *LP,
  467.     int firstEntry, int lastEntry)
  468. /*
  469. This just calls GDSetEntries() or GDDirectSetEntries() to load your ColorSpec
  470. table into the clut of your screen device. It is here simply to provide a
  471. cosmetic match to the call to SetLuminances(), for which loading the clut is
  472. optional. Note: if you prefer, instead of LP you may send just the address of a
  473. ColorSpec table, cast to (luminanceRecord *), since a luminanceRecord begins
  474. with a ColorSpec table.
  475. */
  476. {
  477.     int error;
  478.     
  479.     switch((*device)->gdType){
  480.     case fixedType:
  481.         printf("LoadLuminances: this device has a fixed CLUT.\n");
  482.         return;
  483.     case clutType:
  484.         error=GDSetEntries(device,firstEntry,lastEntry-firstEntry
  485.             ,&LP->table[firstEntry]);
  486.         break;
  487.     case directType:
  488.         error=GDDirectSetEntries(device,firstEntry,lastEntry-firstEntry
  489.             ,&LP->table[firstEntry]);
  490.         break;
  491.     }
  492.     if(error) printf("\007LoadLuminances: error %d\n",error);
  493. }
  494.  
  495. double LToL(luminanceRecord *LP,double L)
  496. {
  497.     if(L > LP->LMax) return LP->LMax;
  498.     if(L < LP->LMin) return LP->LMin;
  499.     return L;
  500. }
  501.  
  502. double VToL(luminanceRecord *LP,double V)
  503. /*
  504. Return the luminance that would result from a nominal voltage.
  505. */
  506. {
  507.     return FixToDouble(_VToL(LP,DoubleToFix(V)));
  508. }
  509.  
  510. FIXED _VToL(luminanceRecord *LP,FIXED _V)
  511. /*
  512. Return the luminance that would result from a nominal voltage.
  513. */
  514. {
  515.     register int i;
  516.     register luminanceTable *LT;
  517.     register FIXED _di;
  518.     double LF,VF;
  519.     
  520.     LT=&LP->L;
  521.     if(LT->exists==luminanceSet){
  522.         _di=FixDiv(_V-LT->_VMin,LT->_dV);
  523.         i=FixToLong(_di);
  524.         _di -= LongToFix(i);
  525.         if(i<0)return LT->_L[0];
  526.         if(i>=LUMINANCES_IN_TABLE-1)return LT->_L[LUMINANCES_IN_TABLE-1];
  527.         return LT->_L[i]+FixMul(_di,LT->_L[i+1]-LT->_L[i]);
  528.     }
  529.     VF=FixToDouble(_V);
  530.     if(LP->powerError < 2.0*LP->polynomialError) LF=VToLPower(LP,VF);
  531.     else LF=VToLPolynomial(LP,VF);
  532.     return DoubleToFix(LF);
  533. }
  534.  
  535. double LToV(luminanceRecord *LP,double L)
  536. {
  537.     return FixToDouble(_LToV(LP,DoubleToFix(L)));
  538. }
  539.  
  540. FIXED _LToV(luminanceRecord *LP,FIXED _L)
  541. /*
  542. New, faster version uses a tabulated version of the gamma function VToL(). This
  543. replaces the old function LToVFormulaic(), formerly called LToV(). A subtlety is
  544. that I tabulate the gamma function _L[]=VToL() rather than its inverse, because
  545. this yields a much lower upper bound on the error of the linearly interpolated
  546. L. This is true because the gamma function is roughly parabolic. Equally spaced
  547. samples of a parabolic function yield equal maximum error of linear
  548. interpolation in all intervals. The cost of tabulating _L[] rather than its
  549. inverse is that we must search the table in order to find the largest _L[] less
  550. than or equal to L. The search time is minimized by using a search procedure
  551. that begins the search at the place found by the last search, since successive
  552. calls to LToV() are likely to request similar luminances.
  553. */
  554. {
  555.     register luminanceTable *LT;
  556.     register int lo,hi,i;
  557.     register FIXED _dL,_dLTemp;
  558.     FIXED _V,_LLo;
  559.     
  560.     LT=&LP->L;
  561.     if(LT->exists != luminanceSet){
  562.         /* get domain of the monotonic section of the gamma function */
  563.         LT->_VMin=DoubleToFix(LToVFormulaic(LP,LP->LMin));
  564.         LT->_VMax=DoubleToFix(LToVFormulaic(LP,LP->LMax));
  565.         LT->_dV=(LT->_VMax-LT->_VMin)/(LUMINANCES_IN_TABLE-1);
  566.         _V=LT->_VMin;
  567.         for(i=0;i<LUMINANCES_IN_TABLE;i++,_V+=LT->_dV) LT->_L[i]=_VToL(LP,_V);
  568.         /* just in case, impose monotonicity, from center on out */
  569.         for(i=LUMINANCES_IN_TABLE/2;i>0;i--)
  570.             if(LT->_L[i-1] > LT->_L[i]) LT->_L[i-1]=LT->_L[i];
  571.         for(i=1+LUMINANCES_IN_TABLE/2;i<LUMINANCES_IN_TABLE;i++)
  572.             if(LT->_L[i-1] > LT->_L[i]) LT->_L[i]=LT->_L[i-1];
  573.         #if FAST_LUMINANCE    /* compute the slope dL/dV */
  574.             /* find shift that maximizes ΔL precision without overflow */
  575.             _dL=0;
  576.             for(i=0;i<LUMINANCES_IN_TABLE-1;i++){
  577.                 _dLTemp=LT->_L[i+1]-LT->_L[i];
  578.                 if(_dL<_dLTemp)_dL=_dLTemp;
  579.             }
  580.             LT->LShift=Log2L(LONG_MAX/_dL);
  581.         #endif                
  582.         LT->exists=luminanceSet;
  583.         LT->latestIndex=-2;    /* invalid, so hunt will start from scratch */
  584.     }
  585.     /* hunt for L in LT->_L[] */
  586.     /* first check at latestIndex, latestIndex±1, ... */
  587.     lo=-1;
  588.     hi=LUMINANCES_IN_TABLE;
  589.     i=LT->latestIndex;
  590.     if(i>lo && i<hi){
  591.         if(_L>LT->_L[i])
  592.             {lo=i;i++;}
  593.         else
  594.             {hi=i;i--;}
  595.     }
  596.     /* simple bisection search, see Numerical Recipes in C, pages 98-99. */
  597.     while(hi-lo>1){
  598.         i=(hi+lo)>>1;
  599.         if(_L>LT->_L[i])lo=i;
  600.         else hi=i;
  601.     }
  602.     LT->latestIndex=lo;
  603.     if(lo<0){
  604.         LT->latestIndex=0;            /* nearest legal index */
  605.         return LT->_VMin;
  606.     }
  607.     if(lo>=LUMINANCES_IN_TABLE-1)return LT->_VMax;
  608.     _LLo=LT->_L[lo];
  609.     _dL=_L-_LLo;
  610.     #if FAST_LUMINANCE
  611.         return LT->_VMin+lo*LT->_dV
  612.             +(_dL<<LT->LShift)/((LT->_L[lo+1]-_LLo<<LT->LShift)/LT->_dV);
  613.     #else
  614.         return LT->_VMin+LT->_dV*(lo+_dL/(LT->_L[lo+1]-_LLo));
  615.     #endif
  616. }
  617.  
  618. double LToVFormulaic(luminanceRecord *LP,double L)
  619. /*
  620. Find a nominal voltage that would give luminance L. Failing that, it returns the
  621. closest possible value. The answer is computed by inverting the formula for the
  622. gamma function.
  623. */
  624. {
  625.     if(LP->powerError < 2.0*LP->polynomialError) return LToVPower(LP,L);
  626.     else return LToVPolynomial(LP,L);
  627. }
  628.  
  629. double VToLPolynomial(luminanceRecord *LP,double V)
  630. /*
  631. Return the luminance that would result from a nominal voltage. Uses a polynomial
  632. equation L = p[0] + p[1]*V^1 +... of any degree. dgp 8/11/89
  633. */
  634. {
  635.     register double L,VV;
  636.     register int i,m;
  637.     
  638.     m=LP->coefficients;
  639.     L=0.0;
  640.     VV=1.0;
  641.     for(i=0;i<m;i++){
  642.         L+=LP->p[i]*VV;
  643.         VV*=V;
  644.     }
  645.     return L;
  646. }
  647.  
  648. double VToLPower(luminanceRecord *LP,double V)
  649. /*
  650. Return the luminance that would result from a nominal voltage.
  651. Uses a rectified power law:
  652.  
  653.     L(V)=power[0]+Rectify(power[1]+power[2]*V)^power[3]
  654.     
  655.     Rectify(x)=x if x≥0
  656.     Rectify(x)=0 if x<0
  657.  
  658. dgp 10/28/89
  659. */
  660. {
  661.     double L;
  662.     
  663.     L=LP->power[1]+LP->power[2]*V;
  664.     if(L>0.0) L=LP->power[0]+pow(L,LP->power[3]);
  665.     else L=LP->power[0];
  666.     return L;
  667. }
  668.  
  669. double LToVPolynomial(luminanceRecord *LP,double L)
  670. /*
  671. Find a nominal voltage that would give luminance L. Failing that, it returns the
  672. closest possible value. Solves a polynomial equation L = p[0] + p[1]*V^1 +... of
  673. any degree. First we find a quick approximate answer by solving the power law
  674. fit, LToVPower(). Then we use Newton's method to home in quickly on the solution
  675. to the polynomial fit. Newton's method is taken from Numerical Recipes in C. For
  676. a polynomial of degree 8 (which is what I recommend) this routine now does about
  677. 3? iterations and takes about 1.1? ms. The error in V is less than TOLERANCE,
  678. i.e. relative to full scale of 256 we get an accuracy of 14.6 bits. You can
  679. increase the TOLERANCE to reduce the number of iterations to make it slightly
  680. faster. A TOLERANCE of 1.0 (for 8-bit accuracy) takes 0.8? ms, not much of a
  681. savings in time, and a large loss in accuracy. Reducing TOLERANCE to 1e-8 yields
  682. an accuracy of nearly 35 bits and takes 1.5? ms, only slightly longer.
  683.  
  684. This routine will be called essentially once per clut entry, so computing a
  685. whole new clut with entries 0..255 will take 256*1.1? ms = 282? ms. Some
  686. experiments use only part of the clut, and will therefore take less time. In
  687. some of the experiments that do use the whole clut it may be desirable to
  688. compute the clut table in advance.
  689. dgp 8/11/89
  690. */
  691. {
  692.     #define TOLERANCE 0.01                    /* desired accuracy of V, see above */
  693.     #define JMAX 20
  694.     double a[MAX_COEFFICIENTS];
  695.     register int i;
  696.     int m,j;
  697.     register double f,df,u,VV,V,dV;
  698.     
  699.     if(L>LP->LMax) return LP->VMax;
  700.     if(L<LP->LMin) return LP->VMin;
  701.  
  702.     m=LP->coefficients;
  703.     if(m>MAX_COEFFICIENTS || m<1)
  704.         PrintfExit("LToVPolynomial: %d coefficients is too many or too few\007\n",m);
  705.     for(i=0;i<m;i++) a[i]=LP->p[i];
  706.     a[0]-=L;
  707.     if(LP->powerError < 0.2*LP->LMax) V=LToVPower(LP,L);    /* very good initial guess */
  708.     else {
  709.         if(LP->quadraticError < 0.2*LP->LMax) V=LToVQuadratic(LP,L);/* fair initial guess */
  710.         else PrintfExit("LToVPolynomial: neither power nor quadratic fit is available\007\n");
  711.     }
  712.     for(j=0;j<JMAX;j++){
  713.         /* evaluate function, i.e. the polynomial, and its derivative */
  714.         f=a[0];
  715.         df=0.0;
  716.         VV=1.0;
  717.         for(i=1;i<m;i++){
  718.             u=a[i]*VV;
  719.             df+=i*u;
  720.             f+=V*u;
  721.             VV*=V;
  722.         }
  723.         dV=f/df;    /* apply Newton's method */
  724.         if(V-dV<LP->VMin) dV=(V-LP->VMin)/2.0;    /* bisect */
  725.         if(V-dV>LP->VMax) dV=(V-LP->VMax)/2.0;    /* bisect */
  726.         V -= dV;
  727.         if(fabs(dV) < TOLERANCE) return V;
  728.     }
  729.     printf("LToVPolynomial(%7.2f): Warning, too many iterations. VToL(%5.1f)=%7.2f\n",L,V,VToL(LP,V));
  730.     return V;
  731.     #undef TOLERANCE
  732.     #undef JMAX
  733. }
  734.  
  735. double LToVPower(luminanceRecord *LP,double L)
  736. /*
  737. Find a nominal voltage that would give luminance L. Failing that, it returns the
  738. closest possible value. Solves the rectified power law:
  739.     L(V)=power[0]+Rectify(power[1]+power[2]*V)^power[3]
  740.     Rectify(x)=x if x≥0
  741.     Rectify(x)=0 if x<0
  742. LToVPower takes about 300 microseconds on a Mac II with 8881 arithmetic. This
  743. routine will be called essentially once per clut entry, so computing a whole new
  744. clut with entries 0..255 will take 256*0.3 ms = 77 ms. Some experiments use only
  745. part of the clut, and will therefore take less time. In experiments that do use
  746. the whole clut it may be desirable to compute the clut table in advance.
  747. dgp 10/28/89
  748. */
  749. {
  750.     double V;
  751.     
  752.     if(L<LP->power[0]) L=LP->power[0];
  753.     V=(pow(L-LP->power[0],1.0/LP->power[3])-LP->power[1])/LP->power[2];
  754.     if(V>LP->VMax) V=LP->VMax;
  755.     return V;
  756. }
  757.  
  758. double LToVQuadratic(luminanceRecord *LP,double L)
  759. /*
  760. This is a quick, approximate version of LToV() that is used by LToVPolynomial
  761. for an initial guess if no power law fit is available. Find nominal voltage that
  762. would give luminance L. Since a quadratic fit to the luminance calibration is
  763. only a fair approximation, this routine serves only to find a quick approximate
  764. answer to the root of a higher order polynomial fit. Solves the quadratic
  765. equation by method recommended in Numerical Recipes in C. My choice of root
  766. assumes that L is an increasing function of V over the operating range of the
  767. display. If concave (i.e. q[2]<0) use the larger root. If convex (i.e. q[2]>0)
  768. use the smaller root. In plain English, if the parabola is right side up then we
  769. take the right hand branch. If the parabola is upside down then we take the left
  770. hand branch. dgp 8/12/89
  771. */
  772. {
  773.     register double a,b,c,q;
  774.     double V;
  775.  
  776.     if(L>LP->LMax) return LP->VMax;
  777.     if(L<LP->LMin) return LP->VMin;
  778.     c=LP->q[0]-L;
  779.     b=LP->q[1];
  780.     a=LP->q[2];
  781.     if(a == 0.0) return -c/b;
  782.     if(b<0.0)
  783.         q=-0.5*(b-sqrt(b*b-4.0*a*c));
  784.     else
  785.         q=-0.5*(b+sqrt(b*b-4.0*a*c));
  786.     if(a*b<0.0)
  787.         V=q/a;
  788.     else
  789.         V=c/q;
  790.     return V;
  791. }
  792.  
  793. double SetLuminanceRange(luminanceRecord *LPtr,double lowLuminance,double highLuminance)
  794. /*
  795. The user will never call this routine directly, as it is called automatically by
  796. SetLuminance and SetLuminances. This routine does the hard work of figuring out
  797. which DACs to fix to optimally represent a given luminance range. The goal is to
  798. fix the coarser DACs and vary the finer DACs to represent the image. The CHANGES
  799. in luminance (i.e. the contrast) will have the accuracy of the coarsest of the
  800. variable DACs. To avoid the overhead of figuring all this out every single time
  801. you call SetLuminance(), the results are stored temporarily in your
  802. luminanceRecord, and are only recomputed if you request a different luminance
  803. range.
  804.  
  805. This function checks and returns immediately if you call it with the same
  806. luminance range that it already has stored in your luminanceRecord.
  807. */
  808. {
  809.     register luminanceRecord *LP=LPtr;
  810.     register int i;
  811.     double dL,dL2,LOffset;
  812.     double VSum,VFixed;
  813.     int dacs,fixed,dac[DACS];
  814.     double gain[DACS],V[DACS];
  815.     double VGoal;
  816.     double lowV,highV;
  817.     double tolerance;
  818.     double VHalfStep;
  819.     double dV;
  820.     unsigned short *RGBArray;
  821.     
  822.     /* -1. If necessary, swap low & high */
  823.     if(lowLuminance > highLuminance){
  824.         dL=lowLuminance;
  825.         lowLuminance=highLuminance;
  826.         highLuminance=dL;
  827.     }
  828.  
  829.     /* 0. Return right away if our work has already been done */
  830.     if(LP->rangeSet == luminanceSet &&
  831.         LP->lowLuminance==lowLuminance &&
  832.         LP->highLuminance==highLuminance)
  833.             return LP->tolerance;
  834.  
  835.     /* 1. sort the DACs by gain, discard any with zero or negative gain.
  836.     Gain is defined as the change in V when a single DAC is increased from 0 to 1.
  837.     Sort the gains so that gain[0]≥gain[1]≥gain[2]...
  838.     Note this order is opposite to that used in Pelli & Zhang (1991). The answers
  839.     are the same, but the notation is different.
  840.     */
  841.     if(LP->rangeSet != luminanceSet){
  842.         /* This only needs to be done once, even if the range changes. */
  843.         gain[0]=LP->r;
  844.         gain[1]=LP->g;
  845.         gain[2]=LP->b;
  846.         #if DACS>3
  847.             #error "Current implementation doesn't support more than 3 DACs."
  848.         #endif
  849.         for(i=0;i<DACS;i++)dac[i]=i;
  850.         Sort(DACS,gain,dac);            /* sort so that gain[0]≥gain[1]≥gain[2] */
  851.         dacs=0;
  852.         for(i=0;i<DACS;i++) if(gain[i]>0.0)dacs++;
  853.         VHalfStep=gain[dacs-1]/2.0;        /* half the smallest step */
  854.         for(i=0;i<COLORS;i++)LP->table[i].value=i;
  855.         LP->dacSize=Log2L(LP->VMax*2-1);    // in bits, rounded up
  856.         if(LP->dacSize!=8)printf(
  857.             "Caution: Your video card seems to have %ld-bit video dacs, and Luminance.c\n"
  858.             "has only been tested with 8-bit dacs.\n",LP->dacSize);
  859.         /* Apple's convention is to provide a 16-bit value to the dac. Luminance.c
  860.         saves time and space by computing values with only as much precision as
  861.         needed for the the video card's dac (i.e. 0 to VMax, or dacSize bits). This
  862.         value must be shifted left by leftShift bits to yield a standard 16-bit value.
  863.         Note that Apple normally recommends scaling the value to fill the entire
  864.         16-bit range, e.g. multiplying an 8-bit value by 0x101, but that doesn't
  865.         seem appropriate here, where we know the dac size, and are more concerned with
  866.         step size that with range.
  867.         */
  868.         LP->leftShift=16-LP->dacSize;
  869.     }
  870.     else {
  871.         for(i=0;i<DACS;i++){
  872.             gain[i]=LP->gain[i];
  873.             dac[i]=LP->dac[i];
  874.         }
  875.         dacs=LP->dacs;
  876.         VHalfStep=LP->VHalfStep;
  877.     }
  878.     for(i=0;i<DACS;i++)V[i]=0.0;    /* initialize fixed dac voltages */
  879.     
  880.     /* 2. Transform lowLuminance and highLuminance to lowV and highV. */
  881.     lowV=LToV(LP,lowLuminance);        /* takes 0.15 ms */
  882.     highV=LToV(LP,highLuminance);    /* takes 0.15 ms */
  883.  
  884.     /* 3. Designate the finest dacs as variable until we have enough to cover
  885.     the range, then designate the rest as fixed. (DACs 0..fixed-1 will be fixed.)
  886.     */
  887.     VGoal=fabs(highV-lowV);
  888.     VSum=0.0;
  889.     fixed=0;
  890.     for(i=dacs-1;i>=0;i--) {
  891.         if(VSum >= VGoal)fixed++;
  892.         VSum += gain[i]*(LP->VMax - LP->VMin);
  893.     }
  894.     /* sum gains of all variable dacs. This is called gVary by Pelli & Zhang (1991) */
  895.     tolerance=0.0;
  896.     for(i=dacs-1;i>=fixed;i--)tolerance+=gain[i];
  897.     /* scale by highest luminance gain in the requested range */
  898.     dL=fabs(VToL(LP,lowV+1.0)-VToL(LP,lowV));
  899.     dL2=fabs(VToL(LP,highV)-VToL(LP,highV-1.0));
  900.     tolerance *= MAX(dL,dL2);
  901.  
  902.     /* 4. Temporarily set the variable DACs to their midpoints. Now set the fixed DACs
  903.     to most accurately represent (lowV+highV)/2.
  904.     */
  905.     VGoal=(lowV+highV)/2.0;
  906.     VSum=0.0;
  907.     VFixed=(LP->VMax + LP->VMin)/2.0;
  908.     for(i=fixed;i<dacs;i++) VSum+=gain[i]*VFixed;    /* The mid-voltage of the variable DACs */
  909.     for(i=0;i<fixed;i++) {
  910.         V[i]=floor((VHalfStep+VGoal-VSum)/gain[i]);    /* the unattenuated voltage of DAC i */
  911.         V[i]=MAX(LP->VMin,MIN(LP->VMax,V[i]));
  912.         VSum += V[i]*gain[i];
  913.     }
  914.     VSum=0.0;
  915.     for(i=0;i<fixed;i++) VSum+=gain[i]*V[i];
  916.     VFixed=VSum;
  917.     
  918.     /* 5. The limited precision of the fixed DACs may result in a 
  919.     small offset between the requested and now available range. The offset will be
  920.     at most half a step of the finest of the fixed DACs. It seems reasonable
  921.     to offset the requested luminance range by up to that small amount to better fit
  922.     the now available range. */
  923.     LOffset=0.0;
  924.     if(fixed>0){
  925.         VSum=VFixed;
  926.         for(i=fixed;i<dacs;i++) VSum += LP->VMin*gain[i];
  927.         dV= VSum - MIN(lowV,highV);
  928.         dV= MIN(dV,gain[fixed-1]/2.0);
  929.         if(dV>0.0) LOffset+=VToL(LP,VSum+dV)-VToL(LP,VSum);
  930.         VSum=VFixed;
  931.         for(i=fixed;i<dacs;i++) VSum += LP->VMax*gain[i];
  932.         dV= MAX(lowV,highV) - VSum;
  933.         dV= MIN(dV,gain[fixed-1]/2.0);
  934.         if(dV>0.0) LOffset-=VToL(LP,VSum+dV)-VToL(LP,VSum);
  935.     }
  936.     /* Now save this information in the luminanceRecord for future use */
  937.     for(i=0;i<DACS;i++){
  938.         LP->dac[i]=dac[i];
  939.         LP->gain[i]=gain[i];
  940.     }
  941.     LP->rangeSet=luminanceSet;
  942.     LP->lowLuminance=lowLuminance;
  943.     LP->highLuminance=highLuminance;
  944.     LP->dacs=dacs;
  945.     LP->VHalfStep=VHalfStep;
  946.     LP->fixed=fixed;
  947.     LP->VFixed=VFixed;
  948.     LP->tolerance=tolerance;
  949.     LP->LOffset=LOffset;
  950.  
  951.     /* Copy into FIXED variables */
  952.     for(i=0;i<DACS;i++)LP->_gain[i]=DoubleToFix(LP->gain[i]);
  953.     LP->_VHalfStep=DoubleToFix(VHalfStep);
  954.     LP->_VFixed=DoubleToFix(VFixed);
  955.     LP->_tolerance=DoubleToFix(tolerance);
  956.     LP->_LOffset=DoubleToFix(LOffset);
  957.  
  958.     /* cache the fixed DAC values, and zero the rest for good measure */
  959.     RGBArray = (unsigned short *) &LP->rgb;    /* treat as an array */
  960.     for(i=0;i<DACS;i++){
  961.         RGBArray[LP->dac[i]] = (short)V[i]<<LP->leftShift;
  962.     }
  963.     return tolerance;
  964. }
  965.  
  966. static void Sort(int n,double arr[],int krr[])
  967. /*
  968. Slightly modified sort routine piksr2() from Numerical Recipes in C. I changed
  969. "float" to "double", and I changed second array to int. I reversed the order, so
  970. largest element would be first. I changed it to use conventional c arrays,
  971. starting at index 0.
  972. */
  973. {
  974.     register int i,j;
  975.     double a;
  976.     int k;
  977.  
  978.     for(j=1;j<n;j++) {
  979.         a=arr[j];
  980.         k=krr[j];
  981.         i=j-1;
  982.         while (i >= 0 && arr[i] < a) {
  983.             arr[i+1]=arr[i];
  984.             krr[i+1]=krr[i];
  985.             i--;
  986.         }
  987.         arr[i+1]=a;
  988.         krr[i+1]=k;
  989.     }
  990. }
  991.